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Ground state properties of multi-orbital Hubbard models are investigated by the auxiliary 
field quantum Monte Carlo method. A Monte Carlo technique generalized to the multi-orbital 
systems is introduced and examined in detail. The algorithm contains non-trivial cases where 
the negative sign problem does not exist. We investigate one-dimcnsional systems with doubly 
degenerate orbitals by this new technique. Properties of the Mott insulating state are quanti- 
tatively clarified as the strongly correlated insulator, where the charge gap amplitude is much 
larger than the spin gap. The insulator-metal transitions driven by the chemical potential 
shows a universality class with the correlation length exponent v ~ 1/2, which is consistent 
with the scaling arguments. Increasing level split between two orbitals drives crossover from 
the Mott insulator with high spin state to the band insulator with low spin state, where the 
spin gap amplitude increases and becomes closer to the charge gap. Experimental relevance of 
our results especially to Haldane materials is discussed. 
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§1. Introduction 

In the last dozen years since the discovery of high Tc materials,0) d and / electron compounds 
have been studied extensively with revived interest for the purpose of clarifying effects of strong 
correlations Ji* In these systems, competition between spin and orbital components plays an impor- 
tant role leading to complex phenomena such as metal- insulator transitions with spin and/or orbital 
orderings and fluctuations. Especially, the relevance of orbital degrees of freedom is a subject of 
recent intensive studies. 

The multi-orbital Hubbard model is one of the realistic and simplified models to investigate the 
interplay between spin and orbital components in the Mott transition. This model is straight- 



1 



forwardly derived from the tight-binding scheme. The expUcit form of the Hamiltonian is given 
by 

H = Ht + 'H£ + + Hu + Ti-j, (1) 



where 



^* = EE*^'^'(4c,v + h.c.) (2) 



ij uu' 



i V 

= -/^ E E "^^^ 
i V 



i u<v' a<(y' 



i v^v' a<a' 

Here the operator c|j^(cj,y) creates(annihilates) an electron in the orbital v = 1, A'^c on the site i = 
1, • • •, A'^s, and riiy = c^^^Ciy defines the number operator. The first three terms are one-body terms: 
TLt is the hopping term, TLf. gives the relative energy level for each orbital and sets the averaged 
chemical potential jj, which controls the electron filling; We take = 0- The last two, Tiu and 

7ij, are two-body interaction terms which are classified as in the conventional way: The former 
term Tiu is often called the on-site Coulomb term which consists of only diagonal interactions. 
The latter 7ij is often called the on-site exchange term which contains the Hund coupling and the 
pair-hopping between orbitals. Note that, in order to keep the rotational symmetry of the two-body 
interaction, the parameters in Eqs. (5) and (6) are not independent of each other originally 
For example, for the twofold Cg orbitals in d electron systems whose wave functions are given by 
il^iia = {x^ - y^) f (r) and Vi2a = ^ (3^^ - r^) / (r), the relation Uu = U22 = Uu + 2Ji2 should 
be satisfied. 

As mentioned above, in these multi-orbital systems, fluctuations of both spin and orbital degrees 
of freedom may compete each other. Various approximations have been examined on this model 
Eq. (|l|); mean-fleld approximations,!' §01' Gutzwiller approximationsJllllll'llllll'0'0' infinite- 
dimensional approache EE) and slave boson approaches. However, these approximations 
in general neglect fluctuations under the competition of various orderings. Although details of 
the mean-fleld approximations in these approaches differ each other, they all have a difficulty in 
describing critical phenomena by nature. A more appropriate treatment which allows studies on 
fiuctuation effects has to be employed in order to discuss critical nature of phase transitions in these 
complex systems. In fact, the criticality may have relevance to many aspects of physical properties 
in c? or / electron compounds. 

In this work, we investigate ground state properties of the multi-orbital Hubbard model by an 
unbiased method, that is, by the auxiliary field quantum Monte Carlo (QMC) method, recently 
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proposed by the authors.E3 Several remarkable aspects which cannot be obtained within one-body 
approximations are elucidated for the one-dimensional systems with doubly degenerate orbitals: 
First, the characteristic properties as the strongly correlated insulator are elucidated quantitatively 
for the Mott insulator at half filling; though the charge and spin excitations both have gaps, the 
charge gap amplitude is much larger than the spin gap. The orbital excitation is gapless, that 
is, the orbital polarization grows continuously with increasing the level split. Critical properties 
of the insulator-metal transition are studied by controlling the chemical potential. Analysis on 
the numerical results based on the hyperscaling hypothesis gives the correlation length exponent 
V = 1/2, which supports the scaling arguments on the Mott transition in one dimension. The 
properties under the self doping between two orbitals are also investigated by controlling the level 
split; the continuous change from the high-spin Mott insulator to the low-spin band insulator is 
elucidated, where the spin gap amplitude increases sensitively depending on the magnitude of the 
self doping and become closer to the charge gap. 

This paper is organized as follows: In the next section, the auxiliary field QMC technique pro- 
posed by the authors is reviewed in detail in a general formalism for arbitrary number of orbital in 
arbitrary dimension. The matrix representation for this formalism and details of QMC updating 
are also described explicitly for the readers who are interested in the actual coding procedure. 
We show that the algorithm does not suffer from the negative sign problem in a non-trivial and 
physically interesting parameter region of the model defined in Eq. (|l]). In Sec. ^, we introduce 
doubly degenerate models to investigate them in the following sections. We also discuss experi- 
mental relevance of this doubly degenerate model. In Sec. the ground state properties of this 
model are investigated by the QMC method in one dimension. Properties of the Mott insulating 
state at half filling is investigated in terms of correlations of charge, spin and orbital degrees of 
freedom. The critical properties of insulator-metal transitions controlled by the chemical potential 
are discussed. For finite level split between two orbitals, effects of self doping are investigated. All 
these QMC results are discussed in Sec. || in detail. Based on the scaling analysis, the critical 
exponents of the insulator-metal transition by filling control are discussed. The crossover between 
the Mott insulator in high spin state and the band insulator in low spin state is examined for the 
level-split control. Finally, Sec. |^ is devoted to summary. 

§2. Method 

Auxiliary field quantum Monte Carlo method is a powerful tool to investigate interacting 
fermions.ii' BiHHliB In this section, we describe a general algorithm for multi-orbital Hubbard 
models defined in Eq. ([l|) for arbitrary number of orbital in arbitrary dimension. This algorithm 
was developed in the previous letter by the authors© We describe it in detail to make this paper 
be self-contained, especially supplementing the matrix representation and the details of updating 
procedure. 
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2.1 Assumption 

For the efficiency of the algorithm described below, we assume Uy^i = [/ in Eq. (5). This allows 
us to factorize the two-body interaction terms Eqs. (5) and (6) into the quadratic forms, as 

^u = ^Y.^n,- ^c? + I (2iVc - 1)E^^ - (7) 

i i 

nj = j:T. ^^L' -EE (8) 

i v<v' i u<u' 

Here Ai^^' = Y^a-^iWu = Y^AAva^iv'^y + Av'a^i^<y)- These factorizations simphfy the Hubbard- 
Stratonovich (HS) transformation as is described below in Eqs. (11) and (12). 

As mentioned in the introduction, the interaction parameters Uyyi and J^yi depend on each other 
to keep the rotational symmetry of the total interaction.i'El* Our assumption Uyyi = U breaks this 
symmetry. In the following, we further treat the values of U and J^u' as independent parameters as 
in many previous theoretical studies. We believe that these may not affect the qualitative features 
of the problem. Note that although our assumptions break the on-site rotational symmetry, the 
symmetry or dimension of elementary excitations might be closely related to the effective exchange 
coupling between different sites which strongly depends on the specific form of the hopping matrix 
tjj' as well as the electron filling or the orbital degeneracy; the following algorithm is applicable to 
any type of t^j' . 

2. 2 Basic formalism 

There exist two schemes in the auxiliary field QMC technique; one is for finite temperatures and 
the other is for the ground state. In the finite temperature calculation, we need to calculate the 
partition function, Z = Tr exp (—/JH) 01311* Compared with this, for the ground state calculation 
which is often called the projection quantum Monte Carlo method,B0^ it is necessary to calculate 
the density matrix, p (/?; (p) = {(p\ exp {—(37i) \4>), where is a trial wave function non-orthogonal 
to the ground state. Below, we describe the basic formalism to calculate the common quantity 
exp {—(3TI) in both schemes. 

First step is the Suzuki- Trotter decomposition into the M = (3 /At slices in the imaginary time 
direction. Then we obtain 

^ (^^-Arm/2^-ArH./2^^Arn,^-Arnu^-ArHj^-ArH./2^-Arm/2Y' ^ q |^^^2^ _ (g) 

Next, in order to derive a one-body expression, we replace the two-body interaction terms 7iu and 
Tij with non-interacting ones by introducing the HS transformation. The general formula of 

the discrete HS transformation we use here is given byiHl 

exp(-ATe/2)= E |exp(isr/iV^/) + 0(Ar4), (10) 

/,s=±l 
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where 7/ = 1 + ^l,r]i = y2{3 — \/6l) and a = At9 > 0. Applying this formula to Eqs. (7) and 
(8), we obtain 



exp(^-Ar|5](n,-iVc)'^ 



n 



exp 



J2 -j-exp{isiir]i^^^/a{{ni- Nc)} 

_hi,sii=±l 



(12) 



nn 

i u<v' 



E 



IVv' ^Vu' L-1 

■'2! '*2i 



exp<^ is2i V'V'^a'' A 



where ai = Ar[//2 > and a^*"' = ArJi,^//2 > 0. 

Of course, the way to decompose the two-body interaction terms, that is, the type of the HS 
transformation is not unique. The advantage of our scheme with Eq. (|l^) lies in the fact that 
it brings the non-trivial parameter regions where the negative sign difficulty does not exist. We 



discuss this point in Sec. 2.5 in detail. 

We apply these decompositions to each Suzuki- Trotter slice in Eq. (^). Then, the one-body 
expression is given as 

M 

exp{-Pn)= ]\J\ [wtcTWeaW^aWUa{h{m),Sl{m))wja{l2{'m),S2irn))WeaWta] (13) 

{/isi/^^'s^"'} ™=1 

with 



Ns 

wta = l[exp{-ATnta/2) 

i=l 

Ns 

i=l V 



W, 



W 



-^+^(2iVc-l)+E^ 



E 



ni, 



wua = 11 -^^^ eyipiisiiiji^^^/a^l^ni^ — 

Ns Jl^ / ^ 

wja= n exp fisai^'V"' V"2'''A!^i''o 



v<v' 



(14) 
(15) 
(16) 
(17) 
(18) 



The product of Eqs. (11) and (12) over all the Suzuki- Trotter slices amounts to the systematic 
error of O(At^) which is negligible because it is higher order than the other systematic errors from 
the Suzuki- Trotter decomposition in Eq. (P). In Appendix A, we discuss the behavior of these 
errors in detail. 
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2.3 Matrix representation 

A matrix representation for the above expressions is helpful in an actual coding. We derive it 
explicitly here in the case of the ground state algorithm for discussions in the following part. The 
following representations can easily be applied to the finite temperature calculation only with slight 
modifications. 

We use a trial wave function |(^) = (g) l^-l-) with 

K) = n E(^<^)ii4 io), (19) 
j=i \i=i ) 

where is the electron number with a spin; iV = iVg x Nq\ and |0) is the vacuum. Then, the 
matrix representation for the density matrix in the canonical calculation is given as 

p (/3; </.) = ({hsxi7's7'}\ 0) ({/isizr'sr'}; /3; 0) , (20) 

{iisi/r'^"'} 

where 

= ^{iisi/r'^"'} '^^^ [*^-^{h.i/r'^r'} ■ ^^^^ 

Here, 

'^{ii.iir'^^'} = 11 11 4 '^''P I "'^^^ ^""^ ^/u(m) 1 (22) 

m=l i=l 

%si;-'3-'}(^2,ri) = n H"^), (23) 

m=mi+l 

where Tj = mjAr (z = 1, 2) and 

6 (m) = e-^^e^e^("^)e^("^)e^e-^^. (24) 



Ei(^fe^. = ^*<^ (25) 



Here T,L,D and are N x N matrices given by 

N 

_ J ^ia "3 
i,j 

{D)Tj = {-m) r]i^.^rn)\/a{5ij (27) 
(^)i] = E i*2^'^' (^) ^^,7'(m) V^^^^- (^-,^^^-.-' + <5.,.5...') . (28) 

In Eqs. (25) ~ (28), the index i contains the site and orbital indices, i = {i, Ui). Sij is the Cronecker's 
delta function. 
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2.^ Updating 

The summation over {lisil2'^ } in Eq. (|20| ) is statistically sampled by the Monte Carlo tech- 
nique. In the ground state algorithm, a change {/, s} = {lisil2^' 82^'} {/', s'} = {I'lS'il"^ s'^ } is 
accepted by the probability 

n.vF.({/',s'}) 



The probability Eq. (|29|), which is called the acceptance ratio, may be rewritten into more 
simple form by using the single-particle Green function.il) For a change from {Ixi (m) , s\i (m)} to 
{Ixi (m) , sxi (m)} (A = 1, 2), the acceptance ratio P is given as 



P 



where / is the N x N unit matrix 

Ci ■■ 

C2 



7«ii(m) 



exp 



Cxl[det [I + AxG. (m)] 



tti^ (sii (m) r/^^^(^) - sii {m)rii^^^^^ 



Aa are sparse N x N matrices with only non-zero elements, 



with 



exp 
exp 



ai isii (m) r?,-^^(^) - su (m) r//^^(„) 



exp 



2 -2^ 



■ / UiUi' ViVi' I \ 



(30) 

(31) 
(32) 

(33) 

(34) 
(35) 



Here, the single-particle Green function at the m-th imaginary time slice is calculated by 

-1 . 



G„ (m) = (m) {^^,B (/3, 0) 



[m 



where 



(m) = e^Me^("^)e-^e-^^^/2^ ((m - 1) Ar, 0) 
*ej (m) = ^^^B {13, {m + 1) Ar) e'^^^/^e^. 



(36) 



(37) 
(38) 



Note that = Cdet [*eje^]. 

If the change of the HS variables is accepted, the updated Green function G for {Ixi {m) , sxi {rri)} 
is calculated from the old one G for {Ixi (m) , sai ("i)} as, 



= (/ + Aa)G,, 



(39) 
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where 

The Nq X A'^c matrix Q is defined as 

V 

To summarize, the actual MC sampling goes as follows: Calculate the single-particle Green 
function for given HS variables. Generate a change of a HS variable and calculate the acceptance 



ratio Eq. (3C) by using the Green function. If the change is accepted, recalculate the Green function 
by Eq. (|39|). Repeat this procedure for all the sites and all the imaginary time slices. This amounts 
to the one MC sweep. 

In general, the weight for each Monte Carlo sample VF^Wj^ can have a negative value, which 
leads to the negative sign problem.ii'0@'ii'i3ili£l@'0) The difficulty caused by the negative 
sign depends on details of the systems; the long-range hopping, the electron density and the di- 
mensionality of the system. In Appendix B, we discuss how the negative-sign samples appear in 



our framework briefly for a specific model which is introduced in Sec. 3.1 



2.5 Cases free from the negative sign problem 

Here we discuss non-trivial conditions on parameters in multi-orbital models which are completely 
free from the negative sign problem.!!^ 

Negative-sign samples do not appear when a system keeps the particle-hole symmetry. There, 
it is easily shown that the weight of the Monte Carlo sample for the up-spin is just the complex 
conjugate of that for the down-spin; W'-fVF^ = |VF^|2 > 0. In multi-orbital cases, the symmetry 
about orbital indices provides us with the following conditions i ~ iv to hold the particle-hole 
symmetry: 

i) The system is at half filling; the electron density satisfies 

n . -^±-^ = No. (42) 

ii) The energy levels split symmetrically around the Fermi energy; 

eu = -ep, (43) 

where u = Nq — v + 1. 

iii) The intra-site exchange couplings satisfy a particular relation; 

Jvv' = Spu'Ju- (44) 

iv) The hopping integrals satisfy a special condition; 

t-' = (-l)l^-^ltg^ (45) 
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where — is the Manhattan distance between site i and j. For these conditions, the particle- hole 
transformation to ensure the positivity of Monte Carlo weights is given by 

c,,t - (-1)^4t (46) 
Ciui Civi ■ (47) 



Here, we illustrate these conditions i ~ iv expressed by Eqs. ( [4 2]) ~ ( [451) and their significance 
by using a heuristic example. We consider here a system with doubly degenerate orbitals in two 
dimensions. An effect of finite e = e\ = —£2 is easily understood in the case of non-interacting 
system with only the nearest neighbor hoppings which are diagonal in orbital indices. Dispersions 
of two orbitals are separated by 2e around the Fermi energy, therefore, a self doping occurs between 
two orbitals. This leads to difference of the Fermi volume and changes of the nesting behavior as 
shown in Fig. 1 (b). Next, we discuss an effect of long-range hopping. We consider here the 
next-nearest neighbor hopping tj^ij^^ = ~'^f{ij)) = ^ 0) where {{ij)) denotes that the sites i and 
j are the next-nearest neighbor pairs. Note that the change in this parameter also controls the 
bandwidth. Fig. 1 (c) shows an effect of this type of long-range hopping. In this case also, a self 
doping between two orbitals occurs. In contrast to the case of a finite e in Fig. 1 (b), it should be 
noted that long-range hoppings can change the shapes of the Fermi surfaces as shown in Fig. 1 (c). 

It should be stressed here that these changes by e or t' do not break the perfect nesting property 
with the nesting vector (vr, vr), because the Fermi surface of the orbital u = 1 coincides with that of 

= 2 by the shift of the momentum (vr, vr). In other words, if the particle-hole symmetry holds, all 
the fc-points on the Fermi surface are nested with those on the Fermi surface of the other orbital. 
This may lead to the Umklapp scattering which may let the system be insulating when we switch 
on the interactions. However, by these changes, the symmetry about the orbital index which exists 
at e = t' = is broken and the self-doping is realized. It is a non-trivial problem how the nature 
of the Mott insulator at e = t' = changes by controlling e or t' . 

Especially, the control of the level split e contains an interesting problem. The system exhibits 
two different states as limiting cases; The Mott insulator at e = and the band insulator at 
£ = 00. The former is a consequence of the strong correlation, that is, it does not appear in the 
non-interacting case. The latter is a trivial insulator which has the same energy gap in both spin 
and orbital channels as the particle-hole excitation gap. It deserves studying how the spin and 



orbital states change with the level split. In Sec. 4.4 and p.2[ we investigate this problem in detail. 



§3. Doubly Degenerate Models 

In the following sections, we concentrate on the doubly degenerate cases of the Hamiltonian Eq. 
(H). By the auxiliary field QMC technique introduced in Sec. the ground state properties of 



doubly degenerate systems are investigated in one dimension. In Sec. 3.1, we define the explicit 
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Fig. 1. Effects of a particular type of the level split or long-range hopping at half filling (See text), (a) Energy contour 
curves for the case with only orbitally-diagonal hopping between nearest neighbor pairs. The arrow represents the 
nesting vector (tt, tt). (b) Energy contour curves for each orbital v = l(left) and 2(right) when the level split 
£ = 0.8t. (c) Energy contour curves for each orbital i/ — l(left) and 2(right) when the next-nearest neighbor 
hopping t' = 0.2t. In all the figures, the Fermi surfaces are denoted by the bold black curves. 



form of the Hamiltonian and give some remarks on basic properties of it. Experimental relevance 



of our model is discussed in Sec. 3.2 



3.1 Hamiltonian and some remarks 

We consider in the following doubly degenerate Hubbard models, given by Eq. (|) with iVc = 2. 
For simplicity, we take into account only nearest neighbor hoppings which are diagonal about orbital 
indices. The level split between orbital u = 1 and 2 is taken to satisfy the condition Eq. (^3|) . 



10 



Moreover, as described in Sec. 2.1, for the convenience of the HS transformation in the auxihary 
field technique, we take uniform values of Ui,^/ irrespective of orbital indices. Then, the explicit 
forms of Hamiltonian are given by 

Wt = -iEEE4.c,- + h.c. (48) 
ne = eJ2{nii-ni2) (49) 

i 

2 

Ti-U = UY^ E E(^~ ^uu'^aa') riiuarii^'a' (50) 
i u<u' = l a<(7' 
2 

i Uy^u' = l u<a' 

where summations on (ij) are over the nearest neighbor sites. In the following, only half-filled cases 
are discussed. Therefore, the particle-hole symmetry holds for all the parameter regions, which 
ensures the absence of the negative sign difficulty in the Monte Carlo calculations as explained in 



Sec. 2.5 



Although this model represents a limited case of the general Hamiltonian (|l]) , it still keeps many 
general aspects and serves our purposes. The hopping diagonal in orbital indices does not cause 
the mixing of two orbitals, while the last term Hj contains off-diagonal elements which mix the 
orbitals v = \ and 2. Note that the second term in Eq. (51) which denotes the pair hopping 
term between two orbitals has been looked over in most of previous works, although this term is 
naturally derived from the tight-binding scheme as mentioned in Sec. ||. 

Our model Eqs. (48) ~ (51) is easily shown to have two insulating states as particular limits: 
(a) When [/, J » t at e = 0, the system is the Mott insulator. The energy gap, which is called the 
Mott gap, becomes 2f7. (b) When e — > oo with finite C/, J and t, since all the electrons fully occupy 
the orbital = 2, the system becomes the band insulator. 

For finite values of U and J, the system is expected to be in the Mott insulating state at e = 
because of the perfect nesting property. In the following, the ground state properties of our model 
are discussed by changing the chemical potential [i and the level split e from this Mott insulator. 
By the former control, we investigate the critical nature of the Mott transition. The latter control. 



from a consideration in the weak correlation regime as detailed in Sec. |2.5| , leads to the self doping 
between two orbitals, which affects spin and/or orbital states in the Mott insulator. We discuss 
how the Mott insulator at e = crossovers into the band insulator at e = c« by this self doping in 
detail. 

3.2 Experimental situation 

In this part, we discuss experimental situations where our model Eqs. (48) ~ (51) may have 
relevance. As mentioned in the previous section, the keywords of our model are two active orbitals 
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and half filling. 

In d electron systems, the fivefold energy levels of the bare d electron is split by the crystal field. 
For example, in the octahedral or tetrahedral surroundings, the fivefold degeneracy splits into the 
twofold eg and the threefold t2g orbitals. For the former octahedral arrangement, the eg levels have 
higher energy than the t2g ones, and for the latter tetrahedral case, upside down. Therefore, the 
on-site situation of our model is realized for the d^ configuration in the octahedral surroundings or 
the configuration in the tetrahedral surroundings. For both situations, the level split between 
the active twofold orbitals is easily induced by the deformation of the surrounding ligands. This 
deformation is considered to be caused by the lattice structure itself or the effect of the external 
pressure. 

There are several materials which consist of the above-mentioned d^ units in the octahedral ligand 
field. One realization is the group of Ni2+ compounds; KNiFg, K2NiF4,i) Ni(C2H8N8)2N02(C104) 
(NENP)i'@) and its related materials or Y2BaNi05.EaE°P Since in the following sections we 
investigate our model in one dimension, we focus here the last two materials, NENP or Y2BaNi05. 
These materials have one-dimensional network of the octahedrons which may be described by our 
model Hamiltonian when we forget the ligands. In these systems, however, the hopping matrix 
is not an orbital-diagonal one like Eq. (48); It may depend on the orbital indices according to 
each lattice structure. Moreover, the level split by the deformation of octahedrons may generally 
take place in a complicated way, not in a symmetric form like Eq. (49). Nevertheless, our simple 
model Eqs. (48) ~ (51) serves to investigate the universal aspects of doubly degenerate systems not 
depending on details of parameters as mentioned in the previous section. In Sec. |5|, our numerical 
results will be examined from this point of view and the relevance on the experimental situations 
will be discussed. 

§4. QMC Results 



In this section, the ground state properties of the doubly degenerate model defined in Sec. |3J. 
are investigated by the auxiliary field QMC technique detailed in Sec. |^. First, numerical details 
in the QMC calculations are mentioned, followed by results of this model. 

4-.1 Details of computations 
The model defined in Sec. |3.l| is investigated in one dimension under the periodic boundary condi- 



tion. The QMC data are obtained from 1, 000 ~ 10, 000 sample averages for the state exp(— ^'H)|(/>), 
where we take /?t = 10 ~ 50 depending on each situation to obtain converged values in the ground 
state. Measurements are divided into five blocks and the statistical error of a quantity is estimated 
by the variance among these five blocks. In order to accelerate the convergence ,0^ we use the 
unrestricted Hartree-Fock ground state as a trial state The QMC calculations are mainly done 
for U = AJ = 2>t. We set Ari = 0.04 for these interaction parameters, which ensures within the 
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statistical errors the convergence to the hmit Ar for all the physical quantities we calculated. 
Detailed discussions of numerical convergence about /3 and At are given in Appendix A. 

4.2 Mott insulating state 

Here we study properties of the Mott insulator at half filling. Low energy excitations in charge, 
spin and orbital degrees of freedom are investigated quantitatively. 

First, we calculate the charge gap amplitude in the Mott insulator at e = 0. The charge gap is 
defined as a change of the total energy when we put an extra particle to the ground state of the 
system with electron number A'e, 

Ac = (iVe + 1) - {Ne) ■ (52) 

In actual QMC calculations, the charge gap is well estimated by the imaginary time dependence of 
a single-particle Green function defined as 

Gyy'cr {rij,T) = {fciya (t) c^^,^ (O) ) , (53) 

where T gives a time ordered product and Ciya (t) = e'^^Qj^o-e"'^^ is in the Heisenberg representa- 
tion. The bracket defines the canonical ensemble average. The Green function in the large r limit 
is governed by the charge gap, which is explicitly shown by the spectral representation as 

lim G^^'^ {rij,T) 

= ^lim |c,..|^^=+^)(v&^^+i|civl^g^=) exp [-r{Em {N, + 1) - E, (N,)}] 

m 

~ «1c...|M/f^+^)(^f=+Vj.vK^) exp [-r{E, (iV, + 1) - E, (N,)}] , 

where I'^m) ^^'^ (-^e) represent eigenvectors and eigenvalues in the A'e electron system respec- 
tively and l^'g) is the ground state. We use the numerically stable algorithm recently proposed by 
Assaad and Imada to calculate the imaginary time dependence of the Green functions.S^ 
Fig. 2 shows typical imaginary-time dependence of a uniform Green function defined as 

G{T) = ^Y.G.uAriuT). (54) 

The charge gap amplitude for each system size is determined by the least square fit for the expo- 
nential tail of G (r). The fitting lines are shown in the figures as the straight lines. The values of 
Ac in the thermodynamic limit are obtained by an extrapolation in the inverse of system sizes as 
shown in Fig. 3 (a). The QMC data are well fitted by the linear function of l/N^- The reason 
for these linear relations is not apparent as it stands. However, this may have a connection with a 
gapless nature in orbital degrees of freedom, which is mentioned below. 

At e = 0, the charge gap as function of interaction strength U/t is summarized in Fig. 3 (b). For 
the comparison, the Hartree-Fock result is also plotted in the figure. Our QMC results show a large 
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reduction from the mean-field results. Our model has both spin and orbital components and there 
is some degeneracy about order parameters. The simple mean-field approach may become crude 
for this type of systems because it tends to favor one ordering among degenerate order parameters 
to open the energy gap and neglect the competition among them. See Appendix C for details of 
the mean-field results. 




Fig. 2. Imaginary time dependence of Green functions at e = 0; (a) for U = — 2t and (b) for U = 4:J — 3t. Note 
that the origin for the y axis has an offset to distinguish the data for each system size in both figures. The straight 
hues are the fit for each data. 



Next, we discuss spin degrees of freedom. In Fig. 4, we show the QMC results of the spin gap As 
calculated as the energy difference between the singlet ground state and the first triplet excitation, 

As = eJn} = ^ + 1;N^ = ^-i)-eJn} = ^;N} = ^). (55) 
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Fig. 3. QMC data for the charge gap. We fix U / J = 4. (a) Extrapolations of the charge gap to the infinite system 
size, (b) The charge gap in the thermodynamic hmit for various values of U/t at e = 0. The gray curve is the 
result of the Hartree-Fock approximation. 

In the figure, we fit the data by As = const. + aN^"^ which was used to estimate the spin gap 
in the S = 1 spin system with nearest-neighbor antiferromagnetic couphng, which is often called 
the Haldane system.© Our data suggest a non-zero value of the spin gap in the thermodynamic 
limit; As/t = 0.12 it 0.05. However, our data are not accurate enough to justify this form of the 
fit, because of large errorbars and finite size limitation. Further investigations are necessary to 
estimate the spin gap As accurately. 

It should be stressed that the magnitude of the spin gap As obtained here is much smaller than 
the charge gap Ac plotted in Fig. 3 (b); As/ Ac 0.1 for [/ = 4J = 4t. This large discrepancy 
between Ac and As is characteristic for the strongly correlated insulator. Note that As is equal 
to Ac in the band insulator. 
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Fig. 4. QMC data for the size dependence of the spin gap at (7 = = At and £ = 0. 



If we take the different values between intra- and inter-orbital Coulomb interactions, that is, 
Uii = U22 > U12, an effective Hamiltonian at e = in the limit of strong correlation is the Haldane 
system. As widely known, this S = 1 spin system has been predicted to have a finite excitation gap 
between the singlet ground state and the triplet excited state, which is called the Haldane gap. Si) 
This conjecture has been favorably supported by both analyticalEl' and numerica 
investigations. Our numerical result obtained above suggests that the Mott insulating state of our 
model at e = may be smoothly connected with the Haldane system, that is, it may be in the 
same universality class in spite of finite values of U and J, although our model with the assumption 
Uii = U22 = U12 = U has an additional degeneracy in the limit of strong correlation. 

Finally, we consider properties in orbital degrees of freedom. In the doubly degenerate system, 
it is useful to describe the two orbital degrees of freedom by the pseudo-spin operator defined as 

^" = ^EE^-'4.c...., (56) 

o" i/u' 

where r is the Pauli matrix and a = x, y or z. In order to discuss the low- lying excitations in this 
channel, we calculate how the pseudo-spin moment responds to its conjugate field, that is, the level 
split e. The pseudo-spin moment is defined as 

-^mom AT I Z_,'^i I' i^'^) 

which signals a magnitude of self doping between two orbitals. As shown in Fig. 5 (a), the pseudo- 
spin moment T^om ^ finite value in the thermodynamic limit when we switch on e. Fig. 5 
(b) summarizes the pseudo-spin moment Tj^om function of e. We find that T^^^ seems to grow 
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Fig. 5. QMC data of the pseudo-spin moment for U = — 3t. (b) shows the values of the pseudo-spin moment 
as function of e, which are obtained from the extrapolations to the thermodynamic limit in (a). 



continuously from zero at e = 0, which strongly suggests a gapless nature in the orbital channel in 
the Mott insulating state. 

4-3 Mott transition by chemical potential control 

We consider here insulator-metal transitions by controlling the chemical potential /j,. Critical 
exponents are investigated from the insulating side by a technique recently proposed and applied 
to single-orbital Hubbard models© In this method, the correlation length exponent u is directly 
determined from the localization length of the single-particle wave function probed by of the 
Green function in the long distance. Here, we explain the prescription to determine the exponent 
u. 

Fig. 6 shows the imaginary time dependence of the Green function for the farthest sites G{rij = 
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Fig. 6. Imaginary time dependence of Green functions for the farthest sites at U = 4J = 3t and £ = 0. 



Ns/2,t) = '^ucjGuuaii'ij = Ns/2,t) for each system size Ns at £ = 0. Prom these data, we obtain 
the frequency dependent Green function for \fi\ < Aq as 

/oo 
dTG(rij,r)e^'^, (58) 

because the density is fixed in the Mott insulating state. It should be noted that the Green function 
satisfies G(rjj,r) = (— 1)1^*^1 G(rjj, —r) by the particle-hole symmetry. In Eq. ([5^), the factor e™ 
makes it difficult to estimate G {rij,uj = fi) accurately, because the statistical error with each QMC 
data will grow exponentially with increasing values of r for tuj > 0. For better numerical estimation, 
we carry out the imaginary time integration in Eq. (58) in the following way@ Using the operator 
Tji which satisfies T^^oiycrTR = CiJ^R^ua, we describe the large r behavior of G {Ns/2,t) as 

G(A^s/2,r)~^(^I/f l4^/2C,..T;Vs/2l^g^=+')(^g^"+VLl<^)exp(-rAc) (59) 

ucr 

in the spectral representation as Eq. (54). Here, our Hamiltonian commutes with the operator 
and we employ the periodic boundary conditions T^^^2 — 1- Therefore, G {Ns/2,t) ~ itG (r) for 
large r, because T;vg/2|^^°) = ±1^'^''). Here G(r) is the uniform Green function defined in Eq. 
(p^). With this property, we estimate the integration in Eq. (|58| ) by substituting the large r 
behavior of G(r), that is, aexp(— rAc), for that of G {Ns/2,t), as 

n+oo r — ro 



(60) 



dT\G (iVs/2, r) le^f" + a e"(^-^c) + / gr(M+Ac) 

-TO \-<J -\-T() J —CO 

Here, we take the threshold value tq such that for r > tq the tail of G {Nq,/2, t) obeys a exp (—rAc) 
within our numerical accuracy. The advantage in Eq. (|60| ) lies in the fact that G (r) generally shows 
better numerical convergence for large r than G {Ns/2,t). 
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Fig. 7 (a) represents A^'s dependence of G{Nq,/2,uj = ^) at e = for various values of the 
chemical potential ji obtained by the above prescription. For each value of /x, G {Ns/2, fi) decays 
exponentially with A'^s due to a finite localization length in the insulating state, that is, G{rij, fj,) ~ 
exp(— |rjj|/^i). We determine the localization length from the decay rates in Fig. 7 (a). Fig. 7 
(b) gives a log-log plot between and | fi — Aq \ , which determines the critical exponent v directly 
through a scaling relation 

6~|/i-Ac|-^ (61) 

We note that the charge gap is determined from Eq. ( |52D independently of this procedure. In this 
case with e = 0, the QMC data may be consistently understood hy u = 1/2. We will discuss this 
result as well as the scaling relation Eq. (^) in Sec. ^.l| . 



4-4 Self doping by level split control 

In this section, we concentrate on the control of the level split e. When we change the level 
split e, the self doping between two orbitals occurs continuously from e = 0, because of the gapless 
nature in pseudo-spin degrees of freedom discussed in Sec. |4.2| . 

Fig. 8 shows the typical behaviors of the momentum distribution function defined as 

"-(^) = ]^E(4.c,..)e'^^-. (62) 

At e = 0, niy=i{k) is equal to n^=2ik) because of the symmetry about the orbital index, and the 
characteristic wave number k* where ni{k) has the largest slope is at 7r/2. When e becomes larger, 
the form of ni{k) changes continuously and the wave number k* approaches zero. Note that the 
relation n2{k) = 1 — ni{TT — k) holds for all the values of e by the particle-hole symmetry. The 
area under the curve of n^{k) corresponds to the density in the orbital u, therefore these behaviors 
of the momentum distribution function clearly represent the self doping from the orbital = 1 to 
u = 2. 



Fig. 9 shows the e dependence of the pseudo-spin moment defined in Eq. (|57|). As shown in Fig. 
9 (a), the pseudo-spin moment T^om grows continuously from zero at e = and approaches 1 for 
large values of e. The behavior of Tj^om large e is also plotted in Fig. 9 (b) as function of e~^. 



This asymptotic behavior of the magnitude of the self doping is discussed in Sec. 5.2 
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Fig. 7. Schemes to calculate the localization length and the correlation length exponent v. See text for details, (a) 
Distance dependence of the Green function for various values of the chemical potential /i. (b) Algebraic divergence 
of the localization length toward the critical point. The gray line corresponds to Ci ~ 1m ^ Ac|~^''^. 



The spin correlation functions for various values of e are shown in Fig. 10. The definition of the 
spin correlation function is 

S{k) = ^Y.^S,S,)e^^''^ (63) 
where 5^ = {S'f + Sf + Sf) /3. The spin operator is defined as 

^r = ^EEC'cLw (64) 

f o-cr' 

with the Pauli matrix r and a = x, y or z. At e = 0, although it has a spin excitation gap as exam- 
ined in Sec. S{k) has a large peak at A; = vr originating from the short-range antiferromagnetic 
correlation in the Mott insulator. When we change the level split e, this structure collapses and 
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k 

Fig. 8. Momentum distribution functions of the orbital 1/ = 1 for the level split control at f/ = 4J = 3t. The gray 
curves are guides to eye. The data are for e/t = 0.0, 0.4, 0.8, 1.0, 1.2 and 1.4 from the top to the bottom. The circles 
and squares correspond to the data for A^s = 14 and 18, respectively. 



the peak is broadened continuously. 

In Fig. 11, we summarize the values of the antiferromagnetic correlation ^(vr) as function of the 
pseudo-spin moment T^om which corresponds to the magnitude of self doping. The inset of Fig. 
11 plots the inverse of S{tt) in the small self-doped region. These behaviors of the spin correlation 
for the self doping are discussed in Sec. ^.2[ 



As mentioned in Sec. 2.5, for a finite value of e, the nesting property between different orbitals 
may play a significant role. Actually, in the mean-field analysis which is discussed in Appendix C, 
the system keeps the insulating nature when we change e from zero, where the energy gap opens 
by the order parameter which corresponds to the scattering process between different spins and 
orbitals. This relevant process is represented by the operator 



^i = \Y.Yl ^u'^a'4uaCiu'a' ■ (65) 



The correlation function of this operator, X(k), is defined as 

X (k) = ^Y.iX,X,)e''^^K (66) 
S ij 

We call this the anomalous correlation function. 

In Fig. 12 (a), we show the behaviors of the anomalous correlation function -'^(A;) with the change 
of the level split e. Note that in the non-interacting case, X{k) takes the constant value 1/2 by 
definition. We plot the values of ^(vr) in Fig. 12 (b). Compared with the monotonic decrease 
of the spin correlation S{tt), X{Tr) is enhanced for small values of e and keeps almost constant 
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Fig. 9. Pseudo-spin moment for the level split control ai U = 4J — 3t. The asymptotic behavior of the self doping 
for large values of e is illustrated in (b). The gray curve is the fit by a/e^ + &/e^. In both figures, the diamonds, 
circles and squares correspond to the data for A'^s ~ 10, 14 and 18, respectively. 



values up to e ~ t. This behavior is discussed in Sec. in comparison with the mean-field results 
obtained in Appendix C. 

§5. Discussion 

In this section, we discuss the numerical results in one dimension obtained in the previous section. 
The Mott transition by the chemical potential control is analyzed based on the scaling theory. The 
behaviors by changing the level split are discussed as the crossover between the Mott insulator 
with high spin state and the band insulator with low spin state. The relevance of our results to 
experimental situations is also remarked. We also compare our results with previous theoretical 
investigations. 
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Fig. 10. Spin correlation lunctions for the level split control at U — AJ — 3t. The gray curves are guides to eye. 
The data are for e/t = 0.0, 0.4, 0.8, 1.0, 1.2 and 1.4 from the top to the bottom. The circles and squares correspond 
to the data for A'^s = 14 and 18, respectively. 
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Fig. 11. The antiferromagnetic spin correlation as function of the pseudo-spin moment at [/ = — 3t. The inset 
shows the behavior of the inverse of ^(Tr) near T^om ~ 0. The circles and squares correspond to the data for 
Ns ~ 14 and 18, respectively. 



5.1 Hyperscaling analysis 

In this part, the QMC data for the insulator-metal transitions controlled by the chemical potential 
obtained in Sec. |4.3| are analyzed based on the hyperscaling hypothesis. The universality class is 
examined by combining scaling relations with the QMC results. 

Note that the hyperscaling analysis assumes the continuous transition. All the numerical data in 
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Fig. 12. (a) Anomalous correlation functions for the level split control at f/ — — 3t. The gray curves are guides 
to eye. The data are for e/t = 0.0, 0.6, 1.0, 1.4 and 2.0 from the bottom to the top at fc = 0. (b) e dependence of 
X{n). In both figures, the diamonds, circles and squares correspond to the data for A'^s = 10, 14 and 18, respectively. 



our investigations do not show any apparent discontinuity expected in first-order transitions. This 
imphes that the transitions are always continuous here, although we can not exclude the possibility 
of a weak first-order transition with small discontinuity. In this section, we assume the continuous 
transitions and employ the hyperscaling hypothesis. 

First, we briefly review scaling arguments for continuous metal-insulator transitionsSillS The 
hyperscaling hypothesis asserts that a single characteristic length scale ^ and a single energy scale 
$7 govern the quantum critical phenomena. The scaling behaviors of these two quantities give two 
basic critical exponents; the correlation length exponent v and the dynamical exponent z, as 

C-IA]-'^ (67) 
^]~e"^~|A|^^ (68) 
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where A is a control parameter for the transition which measures a distance from the critical point. 
With these two exponents v and z, the finite-size scaling for the singular part of the free energy 
density in d dimensions is given by 

/s(A)~A'^('^+^)F(e/L,eV/3), (69) 

where F is a finite-size scaling function, L is a linear dimension of the system and /3 is the inverse 
temperature. Comparing differentiations of /s by L or /? with the Taylor expansion of the free 
energy density for twists of spatial or temporal boundary conditions with dimensional analysis, we 
obtain the scaling behaviors for the singular part of physical quantities; for example, 

charge susceptibility : k ~ A^''^'^^'^^ (70) 
localization length : ^ A^'^. (71) 



The last one Eq. (71) is the same as Eq. (61). 

It is important to notice that the charge susceptibility k can also be deduced directly from the 
definition — 9^//3A^ where A expresses the chemical potential /U. Direct differentiation of Eq. ( |69|) 
gives 

If we compare this relation with Eq. (70), we obtain an important relation 

zv = 1 (73) 
for metal-insulator transitions by the chemical potential. 



In Sec. 4.3, we studied the correlation length exponent v directly from the insulating side. The 
QMC data of the localization length in Fig. 7 (b) suggests v = 1/2. 

The criticality of the filling-control Mott transitions itself in one dimension has been investigated 
based on the scaling theory li There, the critical nature has been discussed depending on whether 
component (spin or orbital) orderings exist or not both in insulators and in metals. The prediction 
of these scaling arguments is that all the Mott transitions by controlling the chemical potential 
in one dimension should be characterized by v = 1/z = 1/2. Our result, v = 1/2 for the doubly 
degenerate Hubbard model also support at least this estimate for v and hence it is plausible that 
it satisfies the hyperscaling with v = 1/ z = 1/2. 

5.2 Effect of self doping 



Here, we discuss the effects of the level split based on the numerical results obtained in Sec. |4.4 . 
In this one-dimensional system, our results indicate the crossover between the Mott insulator with 
high spin state and the band insulator with low spin state with the self doping by the level split. 

At e = 0, the system is in the Mott insulating state. As discussed in Sec. [4.2| , this state has a 
finite spin gap and a gapless excitation in the orbital channel. In this state, two spins in degenerate 
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orbitals couple with each other ferromagnetically by the Hund couphng, that is, they form the 
S = 1 high spin state. These S = 1 spins show the short-range antiferromagnetic correlation which 
leads to the peak structure of the spin correlation function as shown in Fig. 10. 

On the other hand, in the limit e ^ U, J and t, the system should behave as the band insulator. 
Note that for e > 2t, the system becomes an insulator even in the non-interacting case. In this 
limit, two orbitals are split by large e and particles predominantly occupy the one orbital = 2. 
This behavior is clearly seen in Fig. 9 (a). For large values of e, a small density in the orbital 
V = \ may be induced by the 7ij term Eq. (51). By the perturbation in 1/e, the energy gain from 
this mixing is the order of — J^/2e. The fraction of the density in the orbital v = 1 is estimated 
by the e derivative of this gain, therefore its leading term should be 0(e^^). As shown in Fig. 
9 (b), the QMC data of 1 — Tj^om fo'^ large e are well fitted by ae^"^ + be~^ . If we assume the 
rigid-band picture, the energy split between highest-occupied and lowest-unoccupied states is given 
by 2e — W* with the bandwidth W* renormalized by the interactions. Then the e derivative of 
- jV(2e - W*) results in 1 - T^^„^ - ^ + W*/^)^ therefore, the ratio h/a = W* . From the 
fit in Fig. 9 (b), we obtain W* /t = 3.6 it 0.3. This indicates the bandwidth is narrowed by the 
interactions compared with the non-interacting bandwidth 4t in this rigid-band scheme. In the spin 
sector, in this band-insulator limit, the spin moment is close to zero, that is, the low spin state is 
realized. This feature is clearly seen in the flat structure of the spin correlation in Fig. 10. 



Between these two limiting behaviors, our QMC results obtained in Sec. 4.4 suggest the crossover 
from one to the other. Both the pseudo-spin moment and the spin correlation change continuously 
by the level split, as shown in Figs. 9 and 10. In the region e ^ t, as shown in Fig. 12, the 
anomalous correlation deflned in Eq. ( p6| ) has a peak structure at A; = vr. This indicates that the 
Umklapp process between different spins and orbitals is relevant to the insulating nature in this 
region. In the mean-fleld calculations, the importance of this process is also pointed out as discussed 
in Appendix C. Of course, the mean-fleld description is known to break down for one-dimensional 
systems: It predicts the long-range order of this correlation and the phase transition between the 
Mott insulator and the band insulator as shown in Fig. 17. 

When the degenerate orbitals split up by increasing the value of e, our numerical results show 
that the spin correlation is sensitive to a small self doping. As shown in Fig. 11, for [/ = 4J = 3t, 
only 15% self doping is enough to reduce the value of ^(vr) to the half of that at e = 0. For spin 
systems with a finite spin gap, the spin correlation function has a Lorentzian-type peak at A; = vr 
as 

S{k)^[{k-7if + ^fy^'\ (74) 

where is the spin correlation length. This spin correlation length determines the spin gap; 
As ~ ^. Although this form of S{k) is not rigorously applicable to the case for e 7^ 0, it should 
be reasonable to consider that the peak value of S{k) at k = n may be related to the spin correlation 
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length as well as the spin gap at least for e <^ U. The rapid decay of S{tt) against the self doping 
in Fig. 11 suggests the rapid growth of the spin gap from the value in the Mott insulator at e = 0. 
The inset of Fig. 11 indicates that the spin gap increases linearly with increasing the magnitude 
of self doping. 

5.3 Comparison with other experimental and theoretical studies 



In Sec. 5.1, the universality class characterized hy u = 1/z = 1/2 is suggested for the Mott 
transitions controlled by the chemical potential of the doubly degenerate Hubbard model in one 
dimension. This result is considered to be relevant to the metal-insulator transitions when mobile 
holes are introduced into the one-dimensional strongly correlated insulator with a finite spin gap; 
the 5 = 1 Haldane system, the spin-Peierls system or the spin ladder system with even-number 
legs. 

Experimentally, the hole doping to the spin-gap material has been realized in 
(Y,Ca)2BaNi05.i'S) The material Y2BaNi05 is a typical Haldane system with the spin gap 



lOO-ftTjE^EB) where two spins in the e„ orbital in the configuration of Ni^+ construct the one- 



dimensional network as mentioned in Sec. 3^. When we replace with Ca^^ partially, the 
charge doping into this chain is realized. Unfortunately, in this system, the doped holes did not 
show a clean metallic behavior, although the density of states within the Haldane gap is induced by 
them.ii) The temperature dependence of the magnetic susceptibility shows a spin-glass like history 
at low temperature© These experimental results indicate that introduced holes are localized by 
some disorder. For the direct comparison with our result, a material showing the metallic behavior 
with hole doping is desired. Effects of disorder on our theoretical result also remains for further 
study for comparison with the available experimental results. 

Theoretically, several one-dimensional models have been known to have a finite spin excitation 
gap in the Mott insulating state; for example, the t-J-J' model,© the dimerized t-J model,0^ 
the t-J ladder model with even-number legs.0@iB0 These models have been investigated with 
emphasis on the possible persistence of the spin gap in the metallic phase and on a realization of 
the superconductivity near the Mott insulator. From this point of view, the metallic region away 
from half filling in our model provides another interesting example along this line. However, in our 
model, another potential candidate of metals is a ferromagnetic phase due to the double exchange 
mechanism. In this work, we have not investigated this problem mainly because of the negative sign 
difficulty in numerical calculations as discussed in Appendix B. This interesting problem remains 
for further study. 

Next, we give remarks on the experimental relevance of the crossover behavior indicated in Sec. 



5.2 . As discussed in Sec. 3^, there are several one-dimensional systems which are often called 
the Haldane materials. In these systems, the one-dimensional structure itself may induce the 
deformation of the surrounding ligands, which leads to the level split between two orbitals in the Cg 
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state. The magnitude of this level split is difficult to estimate experimentally in general. However, 
it can be the same order as the bandwidth or the Hund coupling. Our results suggest that if this 
happens, the analysis based on the theoretical results for the Heisenberg Hamiltonian may not work 
straightforwardly because of the gapless nature of orbital degrees of freedom as is discussed in Sec. 



4.2. For the finite level split, it is not clear whether the effective spin Hamiltonian justified in the 



strong coupling limit is applicable. Our result implies the necessity to start from the tight-binding 
Hamiltonian as we discussed in this paper at least for quantitative analysis. 



In Sec. 5.2, the linear increase of the spin gap with increase of the self doping is discussed. As 



mentioned in Sec. |3.2| , the self doping may be induced by the change of the ligand field, which is 
caused also by external pressure. Recently, the Haldane material Ni(C2H8N8)2N02(C104) (NENP) 
has been investigated under the external pressure.0^ It shows that the spin gap grows linearly 
with the pressure. Since the pressure does not only simply increase the level split or the self doping 
but also the hopping matrix, that is, the bandwidth, we need a care in analyzing it. However, our 
result on the effect of self doping gives at least one promising explanation for the origin for the 
pressure dependence of the spin gap. 

§6. Summary 

In this work, the multi-orbital Hubbard models are investigated by the unbiased technique, the 
auxiliary field QMC method. The general framework of the numerical method is introduced for 
the multi-orbital models in detail. This technique is applicable for arbitrary number of orbitals 
in arbitrary dimension. We present explicitly that there exist the non-trivial cases where the 
negative sign problem is absent. Using this new technique, the one-dimensional models with doubly 
degenerate orbitals have been investigated in the ground state. Some of basic properties of the 
Mott insulator at half filling are quantitatively clarified. The charge and spin gap amplitudes are 
calculated and the characteristic behavior of them as the strongly correlated insulator is clarified. 
The charge gap can be one order of magnitude larger than the spin gap in our result of the Mott 
insulating phase. The orbital degrees of freedom are shown to have a gapless excitation. By 
controlling the chemical potential, the critical properties of the insulator-metal transition have 
been investigated. We obtain the correlation length exponent i/ = 1/2 which is consistent with 
the general prediction by the scaling arguments. The universality class for this transition has been 
discussed by analyzing our numerical results based on the hyperscaling hypothesis. For the level 
split between two orbitals, the effect of self doping is examined in detail. Our numerical results of 
correlation functions indicate the crossover between the Mott insulator with high spin state and the 
band insulator with low spin state. The significance of the Umklapp scattering between different 
spins and orbitals in this crossover is suggested. With increase of the level split, remarkably different 
spin and charge gap amplitudes in the Mott insulating phase continuously crossover to a similar 
value in the band insulating phase. The antiferromagnetic spin correlation as well as the inverse of 
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the spin gap sensitively decreases with the self doping. Our result of the chemical potential control 
may have relevance to the hole doping of the Haldane materials. Our result of the rapid increase of 
the spin gap by the self doping is consistent with the experimentally observed increase of the spin 
gap in the Haldane materials under pressure. 
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Appendix A 

In this part, numerical convergence in the actual QMC calculations is examined in detail. As 
described in Sec. ^, in the auxiliary field algorithm for the ground state, we should take care of the 
convergence about two quantities; one is the width of the Suzuki- Trotter slice Ar and the other 
is the magnitude of the projection p. We discuss the behavior of these systematic errors for the 



doubly degenerate model defined in Sec. 3.1 



First, we examine the convergence about Ar. As discussed in Sec. ^, the systematic errors 
dominantly come from the Suzuki- Trotter decomposition Eq. (^. The magnitude of the error 
depends on the parameters in the Hamiltonian as well as on the physical quantities which we 
measure. Fig. 13 shows the Ar dependence of the ground state energy per site and orbital, 
Eg = Eg/2Ns, for two sets of parameters, (a) U = 4J = 2t and (b) U = AJ = 3t. The data are 
for 14-site systems with e = 0. Among the quantities which we have calculated, the ground state 
energy has the strongest dependence on Ar. As shown in Fig. 13, the QMC data are well fitted by 

eg = ef -a(Ar)2, (75) 

as expected from the decomposition Eq. (^). Here Sg^* denotes the extrapolated value to the 
limit Ar 0, which is plotted as the diamonds in Fig. 13. In actual QMC calculations, instead 
of obtaining e'^^ by this extrapolation, we take small Ar enough to guarantee the convergence 
within numerical errorbars for each interaction parameters. For example, we employ Art = 0.05 
for U = iJ = 2t and Art = 0.04 for [/ = 4J = 3t. From Fig. 13, these values of Ar are found to 
ensure the numerical convergence to the limit Ar 0. 

Next, the convergence about the projection (3 is discussed. In general, the ground state should be 
projected out from a trial state for large /3 compared with the inverse of the smallest energy scale 
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in the system. Therefore, the values of (3 sufficient to obtain converged data strongly depend on the 
model parameters as well as the trial state. We illustrate this behavior for the change of the level 
split e. Fig. 14 shows the (3 dependence of the pseudo-spin moment defined in Eq. ( |57| ) for various 
values of e. Deviations from the converged values are plotted in the figures. In actual calculations, 
we check the convergence for each parameter and physical quantity carefully, and obtain the results 
by averaging several data for different values of (3 in the converged region. 




Fig. 13. At dependence of the QMC data for the ground state energy of the doubly degenerate Hubbard model. 
The data are for 14-site systems with e = 0; (a) (7 = 4J = 2t and (b) (7 = 4J = 3t. Extrapolated values to At 
by the fit shown as the gray lines are plotted by the diamonds. Filled circles are the data for the values of At 
which we employ in the actual calculations. See text for details. 
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Fig. 14. /3 dependence of the QMC data for the pseudo-spin moment of doubly degenerate Hubbard model for 
various values of the level split. The data are for 14 site systems with U — 4J = 3t. Differences from the converged 
values which we determine are plotted. 



Appendix B 

Here we discuss tlie betiavior of tlie negative sign in our framework introduced in Sec. |2|. Es- 
pecially, we focus two different elements which break the particle-hole symmetry and lead to the 
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negative sign problem; the long-range hopping and the hole or electron doping. 

First, we examine the effect of the long-range hopping. We add the next-nearest neighbor hopping 
proportional to t' to the doubly degenerate model Eqs. (48) ~ (51). Note that this term is 



different from the next-nearest neighbor hopping term discussed in Sec. 2.5; here we consider 



t 



11 



22 



t' > 0. Fig. 15 shows the behavior of the negative sign when we change the value 



of t' at half filling. The negative sign (S) is defined as the real part of the phase of MC weights, 



(5) = Tfe ■ 



(76) 



Next, the effect of the hole doping is discussed. Fig. 16 plots the values of (S) when we introduce 
the extra holes or electrons to the model Eqs. (48) ~ (51) at half filling. We choose the electron 
filling which is closed-shell configuration. As shown in Fig. 16, the negative sign difficulty is far 
severe in this case as compared with the results in Fig. 15. This may be partly because our HS 
transformation Eq. (11) contains the total density as the phase factor; The change of the electron 
filling seems to be crucial to break the phase coherence of each MC samples which is maintained 
at half filling. 



1.0 



<S> 




Fig. 15. /3 dependence of the negative sign for the change of the next-nearest neighbor hopping in 10 site systems 
with U = 4:J — 2t. The dashed lines are guides to eyes. The data are for t' /t — 0.2, 0.4, 0.6, 0.8 and 1.0 from the 
top to the bottom. 



Appendix C 

In this Appendix, by using the mean-field approximation, we discuss overall features of phase 
diagrams for the doubly degenerate Hubbard model defined in Eqs. (48) ~ (51). In Sec. ^, the 
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Fig. 16. P dependence of the negative sign for the hole doping. The filled circles are for TVs = 6 and A^'c = 20. The 
open diamonds are for A^s = 10 and A^c = 12. 



QMC calculations are done for one-dimensional systems. It is difficult to compare the mean-field 
results with the QMC results directly, because strong quantum fluctuations in one dimension lead 
to the breakdown of the mean-field approach. Nevertheless, the mean-field picture is useful to 
understand the basic properties of the model. 

For the chemical potential control, the mean-field phase diagram for the doubly degenerate 
Hubbard model has been investigated by Inagaki and Kubo in three dimensions.!^ They have 
especially concentrated on various spin and orbital states. We discuss here their results near half 
filling which are the issue of our investigation in Sec. |^ and |5[ For e = 0, at half filling, the 
system becomes an insulator with spin antiferromagnetic long-range order and no orbital ordering 
for infinitesimal Coulomb interaction. When we change the chemical potential /x, the insulator- 
metal transition occurs and the metallic state with the antiferromagnetic ordering appears. For 
further change of /x, this spin ordering disappears and the metal without any ordering of spin and 
orbital components is realized. All these features are common to the mean-field phase diagram for 
the single-orbital Hubbard model.El' 

For the change of the level split e, we obtain the Hartree-Fock phase diagrams following the 
method by Inagaki and Kubo.& Here we fix the interaction ratio U/J = 4. Our Hartree-Fock 
calculations are done in one dimension. 

Fig. 17 (a) shows the mean-field phase diagram in the plane of the level split e and the interaction 
strength U . At e = 0, the antiferromagnetic insulator appears for [/ > as mentioned above. In 
this state, the antiferromagnetic order parameter {(^]^^i^Ck+-Kui) has the same value as {c^kyi^Ck+-Ku' i) 
because of the symmetry about the orbital index i/ at e = 0. When we switch on e, the latter 
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{cl^^^Ck+TTu' i) survives and the self doping begins as shown in Fig. 17 (b). There, the insulating 
state persists for e ^ presumably because of the nesting properties between different spins and 
orbitals as discussed in the last part of Sec. ^. The self doping occurs continuously until all the 
electrons occupy the orbital u = 2, that is, the phase transition to the band insulating state. The 
mean-field results for the charge gap Ac which is defined in Eq. ( |5^ ) are plotted in Fig. 17 (c). 
The charge gap has a finite value for all the values of e. However, the origin of the gap is clearly 
different between two phases in Fig. 17 (a). For small values of e, the gap originates from the 
nesting properties and the ordering of {clj^-^Ck+nu' i) , and the value of the gap decreases gradually. 
Compared to this, for large values of e, the gap increases linearly as the function of e, which is 
expected for the band insulating state. 
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